Demographics and reproductive biology of H. schistosus from bycatch on the west coast of India

Analysis for Manuscript

Introduction

This the code for analysis and supplementary figures for the manuscript. The manuscript is currently under review.

#Required libraries

require(tidyverse)
require(lubridate)
require(viridis)

# Importing data

snakes = read.csv("./Data/Sea-snakes_fish-dep_mastersheet_250420.csv")
embryos = read.csv("./Data/Sea snake_embryo data_mastersheet.csv")

snakes <- snakes%>%
  mutate(Year = as.factor(Year))%>%
  filter(Species == "Hydrophis curtus" | Species == "Hydrophis schistosus")%>%
  droplevels()

HS <- snakes%>%
  filter(Species == "Hydrophis schistosus")%>% # filtering out other snakes from bycatch data
  mutate(Year = as.factor(Year))

# function to calculate mcfadden's r2

mcfadden <- function(x){
  
  r2 <- 1 - (x$deviance/x$null.deviance)
  
}

# sampling years
years <- data.frame(Year = factor(c(2016:2019)))

Sample size

The number of snakes (N) sampled from bycatch of trawlers, shore seines and gillnetss.

snakes%>%
  group_by(Species)%>%
  count(name = "N")

Species	N
Hydrophis curtus	282
Hydrophis schistosus	967

The number of snakes sampled in each year:

n.yr <- snakes%>%
  group_by(Species, Year)%>%
  count(name = "N")

n.yr%>%
  spread(Year, N)

Species	2016	2017	2018	2019
Hydrophis curtus	75	38	123	46
Hydrophis schistosus	178	215	521	53

The number of trips (n) and fishing effort (haul-hours) sampled for sea snakes in the current study.

snakes%>%
  filter(Gear.Type != "")%>%
  select(Date, Gear.Type, No..of.Hauls, Average.Haul.Duration..Hours., Tow.duration.hours)%>%
  distinct()%>%
  # Calculating tow duration
  mutate(Tow.duration.hours = case_when(!is.na(Tow.duration.hours) ~ Tow.duration.hours,
                                        is.na(Tow.duration.hours) ~ No..of.Hauls*Average.Haul.Duration..Hours.))%>%
  group_by(Gear.Type)%>%
  summarise(n = n(), # counting number of trips sampled
            haul.hours = sum(Tow.duration.hours, na.rm = T)) # total fishing effort

Gear.Type	n	haul.hours
GillNet	340	535.9633
Rampan	46	190.4500
Trawler	76	285.1600

Demographics of H. schistosus across time

Age structure

# calculating the mean SVL across years

age.yr <- snakes%>%
  group_by(Species, Year)%>%
  summarise(N = n(),
            mean = mean(Snout.to.Vent..cm., na.rm = T))

# SVL cut off for age classes

maturtity <- snakes%>%
  group_by(Species)%>%
  count()%>%
  mutate(juv = 35,
         adult = ifelse(Species == "Hydrophis curtus", 54, 65))

Distribution fo SVL across years

# plotting distribution of SVL across years

snakes%>%
  filter(Snout.to.Vent..cm. > 20)%>%
  ggplot(aes(Year, Snout.to.Vent..cm.))+
  geom_violin(fill = "grey")+
  geom_boxplot(width = 0.1)+
  geom_hline(data = maturtity, aes(yintercept = adult), linetype = "dashed")+
  geom_hline(data = maturtity, aes(yintercept = juv), linetype = "dotted")+
  #stat_summary(fun.data = "mean_sdl", geom = "pointrange", size = 1)+
  geom_label(data = age.yr, aes(Year, 10, label = N))+
  facet_wrap(~Species, ncol = 1, scales = "free_y")+
  labs(y = "Snout to vent length (cm)")+
  theme(strip.text = element_text(face = "italic"))

ggsave(last_plot(), filename = "./Figures/figure1.tiff", height = 6, width = 8)

We did not find any H. curtus neonates in our sampling. Majority of H. curtus were juveniles. H. schistosus populations consisted of mostly adults.

Linear model for change in SVL of each species across years:

library(car)

snakes%>%
  group_by(Species)%>%
  select(Year, Snout.to.Vent..cm.)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(Snout.to.Vent..cm. ~ Year, data = .)),
          sumr = map(mod, broom::tidy),
         stat = map(mod, glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(c(Species:r.squared))

Species	term	estimate	std.error	statistic	p.value	r.squared
Hydrophis schistosus	(Intercept)	85.598303	1.850695	46.2519793	0.0000000	0.0689121
Hydrophis schistosus	Year2017	-13.996389	2.122893	-6.5930743	0.0000000	0.0689121
Hydrophis schistosus	Year2018	-10.158177	1.974836	-5.1438083	0.0000003	0.0689121
Hydrophis schistosus	Year2019	-1.732303	2.818894	-0.6145327	0.5390387	0.0689121
Hydrophis curtus	(Intercept)	65.043404	1.511818	43.0233133	0.0000000	0.1470604
Hydrophis curtus	Year2017	-5.212154	2.405287	-2.1669571	0.0312836	0.1470604
Hydrophis curtus	Year2018	-8.833123	1.825449	-4.8388759	0.0000024	0.1470604
Hydrophis curtus	Year2019	-12.998166	2.225335	-5.8409920	0.0000000	0.1470604

ANOVA to test for effect of years:

snakes%>%
  group_by(Species)%>%
  select(Year, Snout.to.Vent..cm.)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(Snout.to.Vent..cm. ~ Year, data = .)),
          stat = map(mod, car::Anova))%>%
  select(stat)%>%
  unnest()

Species	Sum Sq	Df	F value	Pr(>F)
Hydrophis schistosus	13334.508	3	19.66265	0e+00
Hydrophis schistosus	180165.620	797	NA	NA
Hydrophis curtus	4363.959	3	12.98867	1e-07
Hydrophis curtus	25310.652	226	NA	NA

Plotting fitted values:

snakes%>%
  group_by(Species)%>%
  select(Year, Snout.to.Vent..cm.)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(Snout.to.Vent..cm. ~ Year, data = .)),
         fit = map(mod, ~predict(., years, se.fit = T)),
         fit = map(fit, ~as.data.frame(.)))%>%
  select(fit)%>%
  unnest()%>%
  mutate(Year = factor(c(2016:2019)))%>%
  ggplot(aes(Year, fit))+
  geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
  labs(y = "SVL (cm)")+
  facet_wrap(~Species)+
  theme(strip.text = element_text(face = "italic"))

There was a significant change in SVL distribution of species across years.

Proportion of developmental classes across years

pop.yr <- snakes%>%
  filter(!is.na(Class))%>% # fix 2016 svl
  group_by(Species, Year)%>%
  # Counting number of individuals in each developmental class
  count(Class)%>%
  complete(Class, fill = list(n = 0))%>%
  # Total number sampled in each year
  mutate(N = sum(n))%>%
  group_by(Year, Species, Class)%>%
  # Proportion of each developmental class
  mutate(prop.age = n/N)

# Summary stats

pop.yr%>%
  group_by(Species, Class)%>%
  skimr::skim(prop.age)%>%
  skimr::yank("numeric")%>%
  select(c(Species, Class, mean, sd))

Variable type: numeric

Species	Class	mean	sd
Hydrophis curtus	Adult	0.43	0.31
Hydrophis curtus	Juvenile	0.57	0.31
Hydrophis curtus	Neonate	0.02	NA
Hydrophis schistosus	Adult	0.88	0.09
Hydrophis schistosus	Juvenile	0.14	0.07
Hydrophis schistosus	Neonate	0.03	0.02

Testing change in proportion of developmental classes across years

Fitting a multinomial logistic model to determine change in age classes across years in each species.

require(nnet)
require(broom)

dev.yr <- snakes%>%
  select(Species, Year, Class)%>%
  filter(!is.na(Class))%>% # removing empty class 
  group_by(Species)%>%
  nest()%>%
  mutate(mod = map(data, ~multinom(Class ~ Year, data = .)),
         null = map(data, ~multinom(Class ~ 1, data = .)),
         smry = map(mod, ~broom::tidy(., exonentiate = T)),
         # Calculating mc fadden's r2
         r2 = map2_dbl(mod, null, ~ 1 - .x$value/.y$value))

dev.yr%>%
  select(smry, r2)%>%
  unnest()

Species	y.level	term	estimate	std.error	statistic	p.value	r2
Hydrophis schistosus	Juvenile	(Intercept)	-2.4450606	0.7380571	-3.3128339	0.0009236	0.0455695
Hydrophis schistosus	Juvenile	Year2017	1.1234010	0.7573074	1.4834147	0.1379643	0.0455695
Hydrophis schistosus	Juvenile	Year2018	0.5064641	0.7510767	0.6743174	0.5001096	0.0455695
Hydrophis schistosus	Juvenile	Year2019	-13.3315259	377.0158759	-0.0353606	0.9717922	0.0455695
Hydrophis schistosus	Neonate	(Intercept)	-3.1352560	1.0212895	-3.0698995	0.0021413	0.0455695
Hydrophis schistosus	Neonate	Year2017	-14.5022162	526.2494814	-0.0275577	0.9780149	0.0455695
Hydrophis schistosus	Neonate	Year2018	-0.9347455	1.0900891	-0.8574946	0.3911716	0.0455695
Hydrophis schistosus	Neonate	Year2019	-8.7344420	53.4649902	-0.1633675	0.8702291	0.0455695
Hydrophis curtus	Juvenile	(Intercept)	-1.9386218	1.0661369	-1.8183610	0.0690090	0.0724193
Hydrophis curtus	Juvenile	Year2017	2.7268690	1.1322946	2.4082681	0.0160284	0.0724193
Hydrophis curtus	Juvenile	Year2018	2.5058458	1.0853066	2.3088829	0.0209501	0.0724193
Hydrophis curtus	Juvenile	Year2019	3.5480223	1.1437099	3.1022049	0.0019208	0.0724193
Hydrophis curtus	Neonate	(Intercept)	-11.1507510	99.7723289	-0.1117620	0.9110121	0.0724193
Hydrophis curtus	Neonate	Year2017	-2.2749071	278.6764080	-0.0081633	0.9934867	0.0724193
Hydrophis curtus	Neonate	Year2018	8.2062502	99.7749668	0.0822476	0.9344498	0.0724193
Hydrophis curtus	Neonate	Year2019	-3.4778373	576.2501607	-0.0060353	0.9951846	0.0724193

Plotting fitted estimates of proportions:

dev.yr%>%
  mutate(fit = map(mod, ~data.frame(predict(., years, "probs"))))%>%
  select(fit)%>%
  unnest()%>%
  mutate(Year = factor(c(2016:2019)))%>% # adding years variable
  left_join(n.yr, by = c("Species", "Year"))%>% # getting sample size in each year
  gather(c(Adult, Juvenile, Neonate), key = "Class", value = "p.hat")%>%
  mutate(se = sqrt(p.hat*(1-p.hat)/N))%>% # caluclating standard error
  # Plotting
  ggplot(aes(Year, p.hat, color = Species))+
  geom_pointrange(aes(ymin = p.hat-se, ymax = p.hat+se))+
  labs(y = "Estimated proportion")+
  facet_wrap(~Class)+
  theme(legend.text = element_text(face = "italic"))

Age structure of H. schistosus does not change significantly over a four year period from 2016 to 2018.

The proportion of H. curtus juveniles increases significantly between 2016 and 2018

Change in SVL distribution over months

month.svl <- HS%>%
  filter(Year == "2018", # data not sufficient for other years
         Snout.to.Vent..cm. > 20)%>% # removing erroneous data
  mutate(Month = factor(Month, levels = month.name))%>%
  complete(Month)%>%
  group_by(Month)%>%
  # Calculating mean SVL in each month
  summarise(mean.SVL = mean(Snout.to.Vent..cm., na.rm = T))

# Month of observed births

births <- data.frame(Species = "Hydrophis schistosus", Month = "April")

# plotting distribution of SVL across months

HS%>%
  filter(Year == "2018", # data not sufficient for other years
         Snout.to.Vent..cm. > 20)%>% # removing erroneous data
  mutate(Month = factor(Month, levels = month.name))%>%
  complete(nesting(Species), Month)%>%
  droplevels()%>%
  ggplot(aes(Month, Snout.to.Vent..cm.))+
  geom_violin(fill = "light grey")+
  #geom_point(data = month.svl, aes(x = Month, y = mean.SVL), size = 3)+
  geom_boxplot(width = 0.1)+
  geom_segment(data = births, aes(x = Month, xend = Month, y = 0 , yend = 10), #marking births
               arrow = arrow(length = unit(0.25, "cm"), ends = "first"), size = 1)+
  geom_text(data = births, aes(x = Month, y = 20, label = paste("Observed \n births")))+
  geom_vline(aes(xintercept = "June"), size = 1)+#start of the monsoon ban
  geom_vline(aes(xintercept = "August"), size = 1)+#end of the monsoon ban
  geom_text(aes(x = "July", y = 80, label = "Monsoon ban"))+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  labs(y = "Snout to vent length (cm)")

Proportion of neonates increases in March to May. Birth and neonates were observed around the same time.

Proportion of neonates in each month of 2018

snakes%>%
  filter(Year == "2018", # data not sufficient for other years
         Snout.to.Vent..cm. > 20)%>% # removing erroneous data
  mutate(Month = factor(Month, levels = month.name))%>%
  group_by(Species, Month)%>%
  summarise(prop.neonate = sum(Snout.to.Vent..cm. < 40)/n())%>%
  spread(Month, prop.neonate)

Species	January	February	March	April	May	October	November	December
Hydrophis curtus	0	0	0.50	NA	NA	0	0.0769231	0
Hydrophis schistosus	0	0	0.02	0.05	0.0555556	0	0.0000000	0

Sex ratios

Proportion of females encountered by species

snakes%>%
  group_by(Species)%>%
  filter(Sex == "Male" | Sex == "Female")%>%
  summarise(females = sum(Sex == "Female"),
            N = n(),
            P = females/N)

Species	females	N	P
Hydrophis curtus	102	212	0.4811321
Hydrophis schistosus	304	618	0.4919094

Testing equality in sex ratios

# is proportion of females different from 0.5?

snakes%>%
  filter(Sex == "Male" | Sex == "Female")%>%
  group_by(Species)%>%
  summarise(females = sum(Sex == "Female"),
            N = n())%>%
  group_by(Species)%>%
  nest()%>%
  mutate(test = map(data, ~prop.test(.$females, .$N, p = 0.5)), # Z - test for proportion
         sumr = map(test, broom::tidy))%>%
  select(sumr)%>%
  unnest()%>%
  select(-c(method, alternative))

Species	estimate	statistic	p.value	parameter	conf.low	conf.high
Hydrophis curtus	0.4811321	0.2311321	0.6306857	1	0.4125066	0.5504525
Hydrophis schistosus	0.4919094	0.1310680	0.7173273	1	0.4518628	0.5320580

Sex ratio is not different from 0.5; p = 0.71.

Testing shifts in sex ratios over years

snakes%>%
  filter(Sex == "Male" | Sex == "Female")%>%
  group_by( Species, Year)%>%
  summarise(N = n(),
            females = sum(Sex == "Female"),
            prop.female = females/N)%>%
  ggplot(aes(Year, prop.female, fill = Species))+
  geom_col(width = 0.5, col = "black", position = position_dodge())

Logistic model to test for change in sex ratios of species across years:

snakes%>%
  select(Species, Year, Sex)%>%
  filter(Sex == "Male" | Sex == "Female")%>% # removing empty sex
  mutate(Sex = factor(Sex, levels = c("Male", "Female")))%>%
  group_by(Species)%>%
  nest()%>%
  mutate(mod = map(data, ~glm(Sex ~ Year, family = binomial, data = .)),
         sumry = map(mod, ~tidy(.)),
         r2 = map_dbl(mod, ~mcfadden(.)))%>%
  select(sumry, r2)%>%
  unnest(sumry)

Species	term	estimate	std.error	statistic	p.value	r2
Hydrophis schistosus	(Intercept)	0.5007753	0.2833779	1.7671643	0.0772007	0.0055997
Hydrophis schistosus	Year2017	-0.6120009	0.4379020	-1.3975752	0.1622407	0.0055997
Hydrophis schistosus	Year2018	-0.6052544	0.2977861	-2.0325137	0.0421017	0.0055997
Hydrophis schistosus	Year2019	-0.3404326	0.4010216	-0.8489135	0.3959294	0.0055997
Hydrophis curtus	(Intercept)	1.0185696	0.3235748	3.1478639	0.0016447	0.0593767
Hydrophis curtus	Year2017	-1.3752445	0.5895404	-2.3327403	0.0196618	0.0593767
Hydrophis curtus	Year2018	-1.4885732	0.3812161	-3.9048020	0.0000943	0.0593767
Hydrophis curtus	Year2019	-1.2096248	0.4481189	-2.6993392	0.0069477	0.0593767

Plotting predicted values:

snakes%>%
  select(Species, Year, Sex)%>%
  filter(Sex == "Male" | Sex == "Female")%>% # removing empty sex
  mutate(Sex = factor(Sex, levels = c("Male", "Female")))%>%
  group_by(Species)%>%
  nest()%>%
  mutate(mod = map(data, ~glm(Sex ~ Year, family = "binomial", data = .)),
         fit = map(mod, ~predict(., years, se.fit = T, type = "response")),
         fit = map(fit, ~data.frame(.)))%>%
  select(fit)%>%
  unnest(fit)%>%
  mutate(Year = factor(c(2016:2019)))%>%
  ggplot(aes(Year, fit))+
  geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
  labs(y = "P(Female|Year)")+
  facet_wrap(~Species)+
  theme(strip.text = element_text(face = "italic"))

The proportion of female H. curtus changed significantly across years. H. schistosus sex ratios remained constant across sampled years.

Mortality in bycatch

Overall mortality rate

snakes%>%
  group_by(Species)%>%
  count(Condition.at.encounter..D.A.)%>%
  # Calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)

Species	n	N	prop.dead
Hydrophis curtus	120	282	0.4255319
Hydrophis schistosus	174	967	0.1799380

Testing difference in mortality across species

Logistic model:

snakes%>%
  filter(Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Species, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
  nest()%>%
  mutate(test = map(data, ~glm(Condition ~ Species, family = "binomial", data = .)),
         sumr = map(test, ~broom::tidy(.)),
         r2 = map(test, ~mcfadden(.)))%>%
  select(sumr, r2)%>%
  unnest()

term	estimate	std.error	statistic	p.value	r2
(Intercept)	-0.2814125	0.1209241	-2.327182	0.0199556	0.0502489
SpeciesHydrophis schistosus	-1.2315652	0.1470904	-8.372845	0.0000000	0.0502489

Fitted probability of mortality:

spp = data.frame(Species = c("Hydrophis schistosus", "Hydrophis curtus"))

snakes%>%
  filter(Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Species, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
  nest()%>%
  mutate(mod = map(data, ~glm(Condition ~ Species, family = "binomial", data = .)),
         fit = map(mod, ~predict(., spp, se.fit = T, type = "response")),
         fit = map(fit, ~data.frame(.)))%>%
  select(fit)%>%
  unnest(fit)%>%
  mutate(Species = factor(c("Hydrophis schsitosus", "Hydrophis curtus")))%>%
  ggplot(aes(Species, fit))+
  geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
  labs(y = "P(Dead|Species)")+
  theme(axis.text.x =  element_text(face = "italic"))

H. curtus had as significatly higher mortality rate in bycatch than H. schistosus

Mortality rate by age class

Are juveniles more vulnerable to bycatch mortality than adults?

snakes%>%
  filter(!is.na(Class))%>%
  group_by(Species, Class)%>%
  count(Condition.at.encounter..D.A.)%>%
  #calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)

Species	Class	n	N	prop.dead
Hydrophis curtus	Adult	40	62	0.6451613
Hydrophis curtus	Juvenile	41	125	0.3280000
Hydrophis schistosus	Adult	138	648	0.2129630
Hydrophis schistosus	Juvenile	5	105	0.0476190
Hydrophis schistosus	Neonate	3	8	0.3750000

Testing difference in mortality rates across age classes within species

Logistic model:

snakes%>%
  filter(!is.na(Class),
         Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Species, Class, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
  group_by(Species)%>%
  nest()%>%
  mutate(test = map(data, ~glm(Condition ~ Class, family = "binomial", data = .)),
         sumr = map(test, ~broom::tidy(.)),
         r2 = map(test, ~mcfadden(.)))%>%
  select(sumr, r2)%>%
  unnest()

Species	term	estimate	std.error	statistic	p.value	r2
Hydrophis schistosus	(Intercept)	-1.3012573	0.0960144	-13.5527320	0.0000000	0.0300658
Hydrophis schistosus	ClassJuvenile	-1.6944749	0.4681739	-3.6193277	0.0002954	0.0300658
Hydrophis schistosus	ClassNeonate	0.7904317	0.7365814	1.0731085	0.2832225	0.0300658
Hydrophis curtus	(Intercept)	0.6931472	0.2738613	2.5310156	0.0113733	0.0823342
Hydrophis curtus	ClassJuvenile	-1.3984157	0.3338240	-4.1890813	0.0000280	0.0823342
Hydrophis curtus	ClassNeonate	-16.2592154	1029.1215009	-0.0157991	0.9873946	0.0823342

Fitted probability of mortality across age classes:

dev = data.frame(Class = c("Adult", "Juvenile", "Neonate"))

snakes%>%
  filter(!is.na(Class),
         Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Species, Class, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
  group_by(Species)%>%
  nest()%>%
  mutate(mod = map(data, ~glm(Condition ~ Class, family = "binomial", data = .)),
         fit = map(mod, ~predict(., dev, se.fit = T, type = "response")),
         fit = map(fit, ~data.frame(.)))%>%
  select(fit)%>%
  unnest(fit)%>%
  mutate(Class = factor(c("Adult", "Juvenile", "Neonate")))%>%
  ggplot(aes(Class, fit))+
  geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
  labs(y = "P(Dead)")+
  facet_wrap(~Species)+
  theme(strip.text = element_text(face = "italic"))

H. curtus adults had the highest mortality rate overall. H. schistosus juveniles had the lowest.

Mortality rate by sex

snakes%>%
  filter(Sex != "")%>%
  group_by(Species, Sex)%>%
  count(Condition.at.encounter..D.A.)%>%
  # Calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)

Species	Sex	n	N	prop.dead
Hydrophis curtus	Female	36	102	0.3529412
Hydrophis curtus	Male	52	110	0.4727273
Hydrophis schistosus	Female	62	304	0.2039474
Hydrophis schistosus	Male	72	314	0.2292994

Testing difference in mortality rates across sexes withing species

Logistic model:

snakes%>%
  filter(Sex %in% c("Male", "Female"),
         Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Species, Sex, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
  group_by(Species)%>%
  nest()%>%
  mutate(test = map(data, ~glm(Condition ~ Sex, family = "binomial", data = .)),
         sumr = map(test, ~broom::tidy(.)),
         r2 = map(test, ~mcfadden(.)))%>%
  select(sumr, r2)%>%
  unnest()

Species	term	estimate	std.error	statistic	p.value	r2
Hydrophis schistosus	(Intercept)	-1.3535045	0.1424629	-9.5007493	0.0000000	0.0008565
Hydrophis schistosus	SexMale	0.1453737	0.1957905	0.7424963	0.4577867	0.0008565
Hydrophis curtus	(Intercept)	-0.6061358	0.2071938	-2.9254529	0.0034396	0.0133815
Hydrophis curtus	SexMale	0.5500463	0.2834464	1.9405655	0.0523110	0.0133815

Fitted probability of mortality across sexes:

sex = data.frame(Sex = c("Male", "Female"))

snakes%>%
  filter(Sex %in% c("Male", "Female"),
         Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Species, Sex, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
  group_by(Species)%>%
  nest()%>%
  mutate(mod = map(data, ~glm(Condition ~ Sex, family = "binomial", data = .)),
         fit = map(mod, ~predict(., sex, se.fit = T, type = "response")),
         fit = map(fit, ~data.frame(.)))%>%
  select(fit)%>%
  unnest(fit)%>%
  mutate(Sex = c("Male", "Female"))%>%
  ggplot(aes(Sex, fit))+
  geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
  labs(y = "P(Dead)")+
  facet_wrap(~Species)+
  theme(strip.text = element_text(face = "italic"))

There was no difference in mortality rates across sexes in either species.

Mortality and female reproductive status

Are gravid females more susceptible to bycatch mortality than the rest of our sample?

HS%>%
  filter(Sex == "Female",
         Class == "Adult")%>%
  group_by(Gravid)%>%
  count(Condition.at.encounter..D.A.)%>%
  # Calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)

Gravid	n	N	prop.dead
	41	151	0.2715232
Yes	16	85	0.1882353

Testing thing difference in mortality rate of H. schistosus females by reproductive state

Logistic model:

HS%>%
  filter(Sex == "Female",
         Class == "Adult",
         Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Gravid, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")),
         Gravid = ifelse(Gravid == "", "No", Gravid))%>%
  nest()%>%
  mutate(test = map(data, ~glm(Condition ~ Gravid, family = "binomial", data = .)),
         sumr = map(test, ~broom::tidy(.)),
         r2 = map(test, ~mcfadden(.)))%>%
  select(sumr, r2)%>%
  unnest()

term	estimate	std.error	statistic	p.value	r2
(Intercept)	-0.9685592	0.1834379	-5.280039	0.0000001	0.0087809
GravidYes	-0.4929586	0.3326292	-1.482006	0.1383386	0.0087809

Fitted probability of mortality by reproductive state:

gravid = data.frame(Gravid = c("Yes", "No"))

HS%>%
  filter(Sex == "Female",
         Class == "Adult",
         Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
  select(Gravid, Condition.at.encounter..D.A.)%>%
  mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")),
         Gravid = ifelse(Gravid == "", "No", Gravid))%>%
  nest()%>%
  mutate(mod = map(data, ~glm(Condition ~ Gravid, family = "binomial", data = .)),
         fit = map(mod, ~predict(., gravid, se.fit = T, type = "response")),
         fit = map(fit, ~data.frame(.)))%>%
  select(fit)%>%
  unnest(fit)%>%
  mutate(State = c("Gravid", "Not Gravid"))%>%
  ggplot(aes(State, fit))+
  geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
  labs(y = "P(Dead)")+
  theme(strip.text = element_text(face = "italic"))

Observed reproductive cycle of H. schistosus

Proportion of gravid females in the sample

# percentage of gravid females in sample

HS%>%
  count(Gravid)%>%
  mutate(N = sum(n))%>%
  mutate(prop.gravid = n/N)%>%
  filter(Gravid == "Yes")%>%
  select(-c(Gravid))

n	N	prop.gravid
88	967	0.0910031

Number of gravid females encountered per year

# checking the number of gravid females per year

HS%>%
  group_by(Year)%>%
  filter(Gravid == "Yes")%>%
  count(Gravid)%>%
  spread(Gravid, n)%>%
  rename(`n(Gravid Females)` = Yes)

Year	n(Gravid Females)
2016	1
2018	75
2019	12

Monthly variation in proportion of gravid females

Proper sampling was only conducted in 2018/19 and hence only this period is used for analysis of reproductive cycles.

# calculating the proportion of gravid females per month

gravid <- HS%>%
  mutate(Month = factor(Month, levels = month.name))%>%
  filter(Year == "2018" | Year == "2019")%>% # only for 2018/19
  group_by(Month)%>%
  summarise(N = n(),
            gravid = sum(Gravid == "Yes"),
            prop.gravid = gravid/N)

gravid

Month	N	gravid	prop.gravid
January	102	12	0.1176471
February	134	34	0.2537313
March	129	24	0.1860465
April	54	7	0.1296296
May	36	1	0.0277778
October	16	0	0.0000000
November	69	3	0.0434783
December	34	6	0.1764706

Describing change in the relative proportions of gravid females across months and years of sampling.

# plotting prop gravid per month

gravid%>%
  mutate(Month = factor(Month, levels = month.name))%>%
  complete(Month, fill = list(prop.gravid = 0))%>%
  ggplot(aes(Month, prop.gravid))+
  geom_point(size = 3)+
  geom_path(aes(group = 1), size = 1, linetype = "dashed")+
  geom_segment(aes(x = "April", xend = "April", y = 0 , yend = 0.02), #marking births
               arrow = arrow(length = unit(0.25, "cm"), ends = "first"), size = 1)+
  geom_text(aes(x = "April", y = 0.04, label = paste("Observed \n births")))+
  geom_vline(aes(xintercept = "June"), size = 1)+#start of the monsoon ban
  geom_vline(aes(xintercept = "August"), size = 1)+#end of the monsoon ban
  geom_text(aes(x = "July", y = 0.20, label = "Monsoon ban"))+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  labs(y = "Proportion of gravid females")

Pregnancy from Novemeber to May. Peak in Feb. Observed two live births in April.

Development of embryos and eggs

Sample size

embryos%>%
  summarise(n.mothers = length(unique(Field.Code)),
            n.embryos = length(unique(Embryo.Code)))

n.mothers	n.embryos
29	235

Summary statistics of egg measurements

embryos%>%
  select(Egg.Length..mm.., Egg.Width..mm.., Egg.Weigth..g.., Snout.to.Vent..cm., Embryo.Weight..g.)%>%
  skimr::skim()%>%
  skimr::yank("numeric")%>%
  select(c(skim_variable, mean, sd))

Variable type: numeric

skim_variable	mean	sd
Egg.Length..mm..	39.71	11.72
Egg.Width..mm..	32.23	9.01
Egg.Weigth..g..	14.30	5.45
Snout.to.Vent..cm.	17.72	5.24
Embryo.Weight..g.	4.83	6.28

Sex ratios within clutches

embryos%>%
  group_by(Field.Code)%>%
  filter(Sex != "")%>%
  summarise(prop.female = sum(Sex == "Female")/n())%>%
  skimr::skim(prop.female)%>%
  skimr::yank("numeric")%>%
  select(c(skim_variable, mean, sd))

Variable type: numeric

skim_variable	mean	sd
prop.female	0.53	0.31

Number of feamle embryos

embryos%>%
  filter(Sex != "")%>%
  summarise(females = sum(Sex == "Female"),
            N = n())

females	N
86	166

Testing clutch sex ratio

broom::tidy(prop.test(86, 166, p = 0.5))%>% # Z - test
  select(estimate:conf.high)

estimate	statistic	p.value	parameter	conf.low	conf.high
0.5180723	0.1506024	0.6979603	1	0.4395567	0.5957383

The sex ratio in clutches is not significantly different from 0.5.

Infertility rate

Percentage of eggs with no embryos present:

embryos%>%
  count(Embryo)%>%
  mutate(N = sum(n), rate = n*100/N)%>%
  select(-N)

Embryo	n	rate
Absent	5	2.12766
Present	230	97.87234

Matrotrophy

Do female H. schistosus expend energy in the development of embryos? or Are the eggs formed and yolk only provides nourishment?

Change in yolk weight with embryo weight

embryos%>%
  ggplot(aes(Embryo.Weight..g., Yolk.Weight..g.))+
  geom_point(size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  labs(x = "Embryo weight (g)", y = "Yolk weight (g)")

Linear model to test change in yolk weight with embryo development:

embryos%>%
  select(Yolk.Weight..g., Embryo.Weight..g.)%>%
  nest()%>%
  # Yolk weight is normally distributed
  mutate(mod = map(data, ~lm(Yolk.Weight..g. ~ Embryo.Weight..g., data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:r.squared)

term	estimate	std.error	statistic	p.value	r.squared
(Intercept)	13.0905565	0.2965229	44.14687	0	0.4071222
Embryo.Weight..g.	-0.4526283	0.0450509	-10.04705	0	0.4071222

Yolk weight decreases by 0.45 g for every 1g increase in embryo weight.

Change in egg mass with embryo development

embryos%>%
  ggplot(aes(Embryo.Weight..g., Egg.Weigth..g..))+
  geom_point(aes(col = Yolk.Weight..g.), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  scale_color_viridis(name = "Yolk Weight (g)")+
  labs(x = "Embryo Weight (g)", y = "Total Egg Weight (g)")

Linear model to test relation ship between total egg mass and embryo development:

embryos%>%
  select(Egg.Weigth..g.., Embryo.Weight..g.)%>%
  nest()%>%
  # Total egg mass is normally distributed
  mutate(mod = map(data, ~lm(Egg.Weigth..g.. ~ Embryo.Weight..g., data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:r.squared)

term	estimate	std.error	statistic	p.value	r.squared
(Intercept)	13.5441434	0.4114909	32.91481	0	0.2188866
Embryo.Weight..g.	0.3676922	0.0557912	6.59050	0	0.2188866

Total egg mass increases by 0.36 g for every 1g increase in embryo weight.

Note: Embryo weight is log-normally distributed.

Along with decrease of fat bodies of gravid females (Figure A1) as embryos develop indicates matrotrophic nutrition.

Reproductive effort

Does the amount of enery expended by female H. schistosus reduce with age?

# Cleaning reproductive effort data

re <- embryos%>%
  left_join(snakes, by = c("Date", "Field.Code", "Species"))%>%
  select(Date, Field.Code, Embryo.Code, Species, Egg.Length..mm..:Tail.Length..cm..x,
         Egg.Weigth..g..:Sex.x, Snout.to.Vent..cm..y:Tail.Length..cm..y, Weight..g.,
         -Body.Length..cm..x, - Body.Length..cm..y)%>%
  filter(Species == "Hydrophis schistosus")%>%
  rename(Embryo.SVL = Snout.to.Vent..cm..x,
         Embryo.TL = Tail.Length..cm..x,
         Embryo.Sex = Sex.x,
         Female.SVL = Snout.to.Vent..cm..y,
         Female.TL = Tail.Length..cm..y,
         Mother.Wt = Weight..g.
         )%>%
  # Removing erroneous data point
  filter(Female.SVL > 50)

Increase in clutch size with female age

re_clutch <- re%>%
  select(Field.Code, Egg.Weigth..g.., Mother.Wt, Female.SVL)%>%
  group_by(Field.Code)%>%
  summarise(clutch.size = n(),
            Clutch.wt = sum(Egg.Weigth..g..),
            Mother.wt = last(Mother.Wt),
            Female.SVL = last(Female.SVL))%>%
  mutate(rcm = Clutch.wt/(Mother.wt - Clutch.wt))

re_clutch%>%
  ggplot(aes(Female.SVL, clutch.size))+
  geom_point(size =3)+
  geom_smooth(method = 'glm', method.args = list(family = 'poisson'))+
  labs(x = "Female SVL (cm)", y = "Clutch Size")

Generalised linear model to test change in clutch size with female SVL:

Clutch size is poisson distributed with mean 8.1052632 and variance 10.3216374

re_clutch%>%
  select(Female.SVL, clutch.size)%>%
  nest()%>%
  mutate(mod = map(data, ~glm(clutch.size ~ Female.SVL, data = ., family = "poisson")),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance),
         r.squared =  map_dbl(stat, ~mcfadden(.)))%>%
  select(sumr, stat,r.squared)%>%
  unnest()%>%
  select(c(term:p.value, r.squared))

term	estimate	std.error	statistic	p.value	r.squared
(Intercept)	-0.3764929	0.7174214	-0.5247862	0.5997318	0.560094
Female.SVL	0.0256926	0.0073002	3.5194451	0.0004325	0.560094

The number of eggs borne by females increase by 1.0253151 for every 1 cm increase in female SVL.

Change in overall reproductive effort with age

Beta regression to determine relationship between relative clutch mass and female SVL:

library(betareg)

re_clutch%>%
  select(Female.SVL, rcm)%>%
  nest()%>%
  mutate(mod = map(data, ~betareg(rcm ~ Female.SVL, data = .)),
         sumr = map(mod, ~tidy(., conf.int = T)),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(component:pseudo.r.squared)

component	term	estimate	std.error	statistic	p.value	conf.low	conf.high	pseudo.r.squared
mean	(Intercept)	0.6463332	1.6515322	0.3913537	0.6955358	-2.5906105	3.8832769	0.0497975
mean	Female.SVL	-0.0144713	0.0172664	-0.8381195	0.4019636	-0.0483129	0.0193702	0.0497975
precision	(phi)	12.8451280	4.7017084	2.7320129	0.0062949	3.6299489	22.0603071	0.0497975

re_mod <- betareg(rcm ~ Female.SVL, data = re_clutch)

re_clutch%>%
  mutate(fitted = predict(re_mod, re_clutch),
         q.low = predict(re_mod, re_clutch, "quantile", at = c(0.05)),# lower bound
         q.high = predict(re_mod, re_clutch, "quantile", at = c(0.95)))%>% # upper bound
  ggplot(aes(Female.SVL, rcm))+
  geom_point(aes(col = clutch.size), size = 3)+
  geom_smooth(aes(y = fitted, ymin = q.low, ymax = q.high), stat = "identity", linetype = "dashed")+
  scale_y_continuous(name = "Relative clutch mass")+
  scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
  scale_color_viridis(name = "Clutch size")+
  theme(legend.text = element_text(size = 10))

The overall reproductive effort does not change with female age.

Change in reproductive effort per embryo with female age

Does the effort per embryo change with female age?

# calculating reproductive effort per embruo

re_embryo <- re%>%
  select(Field.Code, Embryo.Code, Egg.Weigth..g..,  Mother.Wt, Female.SVL, Embryo.Sex, Embryo.SVL)%>%
  group_by(Field.Code)%>%
  mutate(clutch.wt = sum(Egg.Weigth..g..), 
         Female.wt = Mother.Wt - clutch.wt)%>%
  group_by(Field.Code, Embryo.Code)%>%
  summarise(inv = Egg.Weigth..g../Female.wt,# investment per embryo
            Female.SVL = last(Female.SVL),
            Embryo.Sex = Embryo.Sex,
            Embryo.SVL = Embryo.SVL)%>%
  ungroup()

Linear model to detemine relationship between female SVL and investment per embryo:

We have controlled for the effect of embryo development.

re_embryo%>%
  select(Female.SVL, inv, Embryo.SVL)%>%
  nest()%>%
  mutate(mod = map(data, ~betareg(inv ~ Female.SVL + Embryo.SVL, data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr,stat)%>%
  unnest()%>%
  select(component:p.value, pseudo.r.squared)

component	term	estimate	std.error	statistic	p.value	pseudo.r.squared
mean	(Intercept)	-2.2798667	0.4894683	-4.657843	0.0000032	0.5627972
mean	Female.SVL	-0.0161182	0.0050085	-3.218169	0.0012901	0.5627972
mean	Embryo.SVL	0.0491682	0.0077622	6.334311	0.0000000	0.5627972
precision	(phi)	233.3499725	43.5079611	5.363386	0.0000001	0.5627972

Plotting relatonship between residuals and female SVL:

# using residuals to control for embryo development

emsvlinv <- betareg(inv ~ Embryo.SVL, data = re_embryo)

# Plotting

re_embryo%>%
  modelr::add_residuals(emsvlinv)%>%
  ggplot(aes(Female.SVL, resid))+
  geom_point(aes(col = Embryo.SVL), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  scale_y_continuous(name = "Relative egg mass (residuals)")+
  scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
  scale_color_viridis(name = "Embryo SVL (cm)")

The relative egg mass (controlled for embryo development) reduces with female SVL.

Difference in reproductive effort for male and female embryos

Does female investment differ by sex of the embryo?

re_embryo%>%
  filter(Embryo.Sex == "Male" | Embryo.Sex == "Female")%>%
  modelr::add_residuals(emsvlinv)%>%
  ggplot(aes(Female.SVL, resid, shape = Embryo.Sex))+
  geom_point(aes(col = Embryo.SVL), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  scale_y_continuous(name = "Relative egg mass (residuals)")+
  scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
  scale_color_viridis(name = "Embryo SVL (cm)")+
  scale_shape_discrete(name = "Embryo Sex")

Linear model to test for the difference in female investment by embryo sex:

We controlled for the effect of embryo development.

re_embryo%>%
  filter(Embryo.Sex == "Male" | Embryo.Sex == "Female")%>%
  select(Embryo.Sex, inv, Embryo.SVL)%>%
  nest()%>%
  mutate(mod = map(data, ~betareg(inv ~ Embryo.Sex + Embryo.SVL, data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(component:p.value, pseudo.r.squared)

component	term	estimate	std.error	statistic	p.value	pseudo.r.squared
mean	(Intercept)	-3.7036341	0.1763377	-21.003076	0.0000000	0.4781387
mean	Embryo.SexMale	-0.1376065	0.0839749	-1.638663	0.1012835	0.4781387
mean	Embryo.SVL	0.0481588	0.0083489	5.768291	0.0000000	0.4781387
precision	(phi)	205.9146533	38.4135707	5.360466	0.0000001	0.4781387

Female investment does not differ significantly with embryo sex.